arXiv:1508.07344v2 [physics.comp-ph] 23 Nov 2015 


A Spectral Canonical Electrostatic Algorithm 


Stephen D. Webb^ 

^RadiaSoft, LLC, 1348 Redwood Ave., Boulder CO 80304 
E-mail: swebb@radiasoft.net 

Abstract. Studying single-particle dynamics over many periods of oscillations is a 
well-understood problem solved using symplectic integration. Such integration schemes 
derive their update sequence from an approximate Hamiltonian, guaranteeing that the 
geometric structure of the underlying problem is preserved. Simulating a self-consistent 
system over many oscillations can introduce numerical artifacts such as grid heating. 
This unphysical heating stems from using non-symplectic methods on Hamiltonian 
systems. With this guidance, we derive an electrostatic algorithm using a discrete 
form of Hamilton’s Principle. The resulting algorithm, a gridless spectral electrostatic 
macroparticle model, does not exhibit the unphysical heating typical of most particle- 
in-cell methods. We present results of this using a two-body problem as an example 
of the algorithm’s energy- and momentum-conserving properties. 
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1. Introduction 

The numerical simulation of plasmas using macroparticle model^ [H [2] is a common 
tool for applications ranging from plasma processing to particle accelerators, fusion 
to laser-plasma interactions. The traditional approach to such algorithms starts with 
discretizing the Lorentz force law in time, assuming the helds are given everywhere in 
space. The Poisson or Maxwell equations are then discretized in space and time, and 
updated synchronous to when the Lorentz force integrator requires them. The helds 
used for the Lorentz force law are interpolated from the discrete spatial points at the 
proper time. Depositing the source terms is then chosen to carefully avoid self-forces 
on a single particle and conserve the local charge density. The update sequence is also 
frequently chosen to conserve momentum, energy, or both. 

This process is prone to unphysical heating due to various instabilities having to 
do with the grid, the time discretization, or both. They also impose conservation laws 
which do not arise naturally from the discretization. Essentially, the problem is that 
when discretizing the equations of motion there are too many degrees of freedom left to 
the algorithm designer which may lead to unphysical behavior. 

This can be remedied by deriving the algorithms from a Hamiltonian least-action 
principle. Such multisymplectic algorithms, described most recently by Shadwick, 
Stamm, and Evstatiev respect the fact that the dynamics of both particles and 

helds [5l [6] have a symplectic geometry. The benehts of symplectic integration for 
single-particle mechanics are well-known in the accelerator physics community (see, for 
example, [7j for an historical overview of the subject and its applications). The idea 
of multisymplectic particle-in-cell algorithms has generated considerable interest, and 
a number of one-dimensional examples on a discrete spatial grid have recently been 
published [H [9l [101 E] 

The chief beneht of multisymplectic algorithms is that their solutions must be 
a solution for some Lagrangian which approximates the continuous, “real world” 
Lagrangian. The numerical solutions are exact solutions to an approximate Lagrangian. 
This contributes to long-time stability over many oscillations of the system. In the 
case of plasmas, such algorithms promise to make simulations extremely reliable over 
hundreds or even thousands of plasma periods. 

By making the discretizing approximations in the action integral, and then deriving 
the resulting equations of motion, we remove many of the undetermined components 
that can appear when approximations are made to the equations of motion directly. For 
example, charge deposition and force interpolation arise from a common term in the 
Lagrangian, and thus any approximations made to the particles and helds will relate 

f We distinguish between a macroparticle model - which represents large numbers of point particles 
with a single particle of finite spatial extent - and a particle-in-cell model - which places those finite 
spatial extent particles into a spatial grid and carries out charge deposition by interpolating to the 
grid. Thus, particle-in-cell models are a subset of macroparticle models. In spectral algorithms, both 
options are viable, although in this paper we consider only gridless methods. 
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these two constructively. The textbook method is to insist on a deposition/interpolation 
scheme which does not produce self-forces on the particles. This may or may not be 
consistent with the underlying physics which produces the equations of motion. 

What is important to understand here is that the familiar physical principles ~ 
conservation of momentum, conservation of energy, conservation of angular momentum, 
and on - are encoded not in the equations of motion but in the symmetries of the 
action. They appear in the equations of motion, but arise because of symmetries under 
translation in space and time, rotations, and on, per Noether’s Theorem. If we want our 
algorithms to contain these conservation laws consistently, then our algorithms should 
come from some action principle which has the appropriate symmetries. 

In this paper we introduce an electrostatic canonical gridless macroparticle 
algorithm using a spectral decomposition of the electrostatic helds. We express the 
Lagrangian in terms of the electrostatic potential and particle phase space density 
(Section]^, then discretize the Lagrangian into Fourier modes for the helds and macro¬ 
particles for the phase space density (Section |^. In Section the resulting action 
integral is discretized into a second order Riemann sum, and we use the discrete Euler- 
Lagrange equations to arrive at a second order in time spectral algorithm. The linear 
numerical dispersion relation is derived in Section We demonstrate the momentum 
and energy conservation properties for a two-dimensional two-body problem in Section]^ 
and for a thermal plasma in Section Finally, we discuss some optimization strategies 
for implementation in Section 

2. The electrostatic Lagrangian 

Single particle symplectic algorithms can be derived either from a sequence of canonical 
transformat ions na, splitting methods or so-called “kick codes” [13], or from a discrete 
action minimization m and references therein]. Regardless of the formalism used 
to derive the symplectic algorithm, the numerical solutions are an exact solution to 
the equations of motion derived from an approximate action. This is a much tighter 
constraint than schemes which start with the equations of motion, where any arbitrary 
discretization scheme likely fails to respect the geometric properties required by a system 
described by a Hamiltonian least-action principle. This is what gives them their long¬ 
term stability. 

In this paper, we consider the Lagrangian treatment of the least-action principle. 
This contains the same geometric structures as a Hamiltonian approach, but will be more 
convenient for future work with electromagnetics, where the canonical momentum’s 
relation to the vector potential makes explicit integration schemes much more difficult. 

We begin with the electrostatic Lagrangian for a system of particles and helds, due 
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to Lo'ftjl] [15] 
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/(xo,vo) + ... 


( 1 ) 


An auxiliary requirement is that the phase space volume is conserved, i.e. /(xq, vq, 0) = 
which is a statement of the Vlasov equation. The continuum equations of 
motion come from minimizing the action S = f dtC with respect to x and (p, and the 
resulting equations are the familiar Newton’s Second Law for the particle characteristics, 
and the Poisson equation: 



(2a) 


(26) 


To derive our discrete equations of motion from an action, we must hrst decompose 
C into a discrete set of coefficients the basis functions of our phase space density and 
scalar potential, and then approximate S as some discrete approximation of the action 
integral. The former is a mode decomposition, while the later is an approximation of 
the action integral using Riemann sums. Once this is done, we will use the discrete 
Euler-Lagrange equations from [IT] to get the discrete equations of motion. 

3. Spatial discretization of the Lagrangian 

Two helds must be discretized - the phase space density and the scalar potential. 
Discretizing tp will specify the Poisson solving algorithm, while discretizing the phase 
space density will give us a dehnition of macroparticles. 

We discretize the scalar potential in a discrete Fourier-type basis, allowing for shape 
functions in Fourier space: 



( 3 ) 


where a is a collective index for the discrete modes and T is a generic shape function 
for the Fourier modes. 

The choice of spectral shape function can be made to minimize ringing. For 
example, a delta function will yield the unattenuated cos{kax) oscillation. By using 
tent functions separated by Ak, on the other hand, it is possible to accurately represent 
the Fourier transform of the helds in a piecewise linear fashion. This puts an envelope 
on the oscillations, and damps out ringing. Other options are also possible in this 
formalism, although we must consider their specihc properties, as we do in § [^ below. 
In this way, the delta function option may prove the fastest, although it has some issues 
with ringing. An example of the tent function shape in real space is given in hgure Q. 


§ We use C.G.S. units throughout this paper. 
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Figure 1: The real and imaginary parts in the x space for a tent shape function — fco-)- 
The envelope drops off as {x Ak)~‘^ for the tent function. 


To discretize the phase space density, we take shape functions in real space for the 
positions, and use delta functions in velocity: 


/(x,v,f) = ^M;jA(x-Xj(f))5(v-Vj(f)) (4) 

j 

where Xj and \j are the coordinates and velocities for the macroparticle, which has 
a weight Wj. Here, we use v as shorthand for c^tx. We can understand macroparticles in 
this context as a finite element Lagrangian picture of the flow of the phase space fluid. 

Inserting these expressions into the action, eqn. Q, yields a discrete Lagrangian 
whose independent dynamical variables are the Fourier mode coefficients of the scalar 
potential and the coordinates and velocities of the macroparticles: 


■ID — 



j dkA*(k)e**^'"^T(k-k.)(p„ 


+ ^ J k ■ kT(k - k^)T(-k - k^/). (5) 

(7 a' 

We have used the orthogonality of the Fourier modes when integrating over x in the 
scalar potential energy term, and the particle shape function Fourier transform, A(k), 
is dehned the same as the Fourier transform for tp above. Specifically, this requires that 

A(x) = f dke'^'"A(k) ^ A(k) = j dxe-'‘^'"A(x) (6) 

from the Dirac delta function identity (5(x) = (27r)^ J (ike*’^ ’^. For convenience we will 
dehne the /c-squared matrix as 


/Co-cr' = 


dk k ■ kT(k - k,)T(-k - k,0- 


(7) 
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If we consider the simple discrete Fourier transform, with d'(k — kg-) = (5(k — ko-), 
this becomes an uncoupled Poisson equation, where modes are coupled only to their 
complex conjugates. This is the functional form of the continuum Lagrangian, which 
is invariant under translations and therefore conserves momentum. For a more general 
shape function, /C has off-diagonal elements that couple different modes. This matrix 
must be inverted once at the beginning of the simulation to compute the helds. It also 
requires evaluating various convolution integrals between the particle and held shapes, 
which may prove computationally prohibitive if an analytical expression is not available. 


4. Discretizing the action and the discrete equations of motion 


The continuous-time action integral is now approximated by a discrete Lagrangian which 
is still a function of continuous time: 
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dth 


£)[Xj, Vj, (Po-J 


( 8 ) 


To hnd the discrete equations of motion, we must approximate the integral for S using 
discrete steps in time, using Riemann sums or higher order approximations. The most 
direct approach is to use composition methods to make a self-adjoint approximation of 
the action (see §2.4 and §2.5 of H)- 

A hrst order approximation of the action is the Riemann sum 


= KLd 


(n+l) _ (n) 

yn+1) , 

’ h ^ 


A, 


(n+l) 


( 9 ) 


for a step size h. This is the action for the step, and the total action is approximated 
by5«E„st”> -|- 0{h‘^). This is a hrst order approximation. To reach second order 
accuracy in time, we use a composite scheme of a step with its adjoint. The adjoint of 
the discrete action is given by 

h) = -Sfl(x("), -h) (10) 


and it is straightforward to show that a self-adjoint action is of an even order of accuracy 
- thus by being self-adjoint we are assured at least second order accuracy in time. 

By taking the composite action for the step 


Sd = ^ ^mwi 




^2 








(n+ 1 / 2 ) 


ptcl. move 

dk AYk)T(k-k, 




(n + l/2) 


interpolation/deposition 


Svr 




,if 


(n+1/2) 


( 11 ) 


Poisson equation 
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and minimizing the action nsing the discrete Enler-Lagrange equations (see §2.5 of iia: 
we get the discrete equations of motion 


^(n+l) 


2x 


(n+1/2) 


m- 


(n) 


Y 2 


Ek 


dk kA*(k)\l/(k — 




^(n+ 3 / 2 ) _ ^(n+ 1 ) 


= Svre Wj 


(n+l/2) 

dkA*(k)e*‘^E ^(k-k. 


j 


An+l) _ ^(n+^/2) 


(12 a) 

(126) 

(12c) 


h /2 hj2 

where /Co-o-' is the generalization of the usual k^ in the Poisson equation. The first 
equation is simply mx = gE, while the second is the Poisson equation for computing 
the potential in the middle of the step. The third equation tells us how the boundaries 
between steps are related. In our case it is the continuity of the “velocity”, although 
different discretizations may have different boundary terms. 

Note that there is no velocity coordinate - as discussed by Marsden and West [H] 
there is no tangent space in the discrete Lagrangian picture. We may, however, introduce 
the velocity as a convenient intermediate computational variable. If we define 

x(’i+V2) _ -An) 

(13) 


(n) _ 


' = 


Y 2 


then the above equations become 


(n+1/2) _ (n) ^An) 

W “W 2 


/Co-o-'^^cr' = Svre ^ tCj J dkA* (k)e**^'’‘^ ^' ^(k — k, 

^ y dk kA*(k)T(k - 


(71+1/2) _ (71) _ 

j j 


J 

he 

—i 

m 




(14a) 

(146) 

(14c) 


(71+1) _ (n+V2) I ^ (n+V2) 

A ~ A ^2 ^ 

( 71 + 1 ) _ (71+1/2) 

j j 

which is the split-step move and accelerate familiar from most particle-based algorithms. 
We have avoided using a grid to deposit charge to preserve the translational invariance. 
This means that every particle has an exponential, that must be calculated. 

We have thus far left T almost completely general. To understand conditions on 
momentum conservation, we derive the total force between two particles. First, let us 
assume that we have chosen our discrete k^. such that there exists a k_o- = —kg.. This 
is a reasonable choice if we want to make the shape functions A and T even functions 
and still return real values for </5(x). 







If we place a single particle at position x^, it will create the potential 




STreWj f — ko 

a' 


where 


/.(x) = / rfkA*(k)^(k-k.)e* 


The force on particle i due to particle j is then given by 

Apj^e = STihe^weWjV^i, ^ /^/(xj)lC“J,/^(x£) (17) 

aa' 

If we dehne 

F(k,k') = A*(k)A*(k')f(k + k')e*(‘^'^^+*^''^^) (18) 

then we can then say that the total change of the momentum for this two-body pair is 

/ OO 

dkdk' F(k, k') T(k - k.)/C;J,T(k' - k^O 

/ OO 

dkdk' F(k, k') T(k - k_,)/C:^_^,T(k' - k_,0 

/ OO 

dkdk' F(k, k') ^ T(k + k^)/C;J,T(k' + k^O (19) 

where we have used to symmetry of )C-a-a' = ^a,a' and the fact that a and a' are 
dummy indices. We can then flip the sign on k and k' and exploit the fact that T is an 
even function to get that 

/ OO 

dkdk' F(-k, -k') T(k - k.)/C;J,T(k' - k^O 

•OO ^ 

(7,<7 

= (Ap,^, + Ap,^,)* (20) 

This proves that the total change in momentum is a real number if we choose even 
functions for our particle and held shapes, and have negative pairs of Fourier modes. 
However, this will not generally be zero, and we can have some change in the total 
momentum. This is because the resulting discrete Lagrangian is not, in general, 
unchanged by translations, so it does not have the associated Noether conserved 
quantity. 

For the special case of a delta function, the Lagrangian is unchanged under a 
translation, and we can show immediately that momentum is conserved to machine 
precision. Here, = k^. • ko-'(5o--o-'. There are no self-forces 


Apj_,.j = Airhe^wi 
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since every is paired with a —ko- = k_a-, a requirement on the expansion to keep </3(x) 
real. Because there are no self-forces, we can write the force on particle £ from particle 
j as 

Api^j = iTihe^weWji ^ = -Apj^i. ( 22 ) 

■^cr 

cr 

We must therefore balance our desire for preventing the ringing from using discrete 
Fourier modes with how much our problem requires exact or nearly-exact momentum 
conservation. It remains to see whether such systems will grow monotonically in 
momentum, or whether the total momentum will simply oscillate around the initial 
total momentum, and how these depend on our shape functions. For now, we will focus 
on the discrete modes with delta-function shape, as they are the least computationally 
demanding due to the lack of any numerical integrals in the update sequence. 


5. Numerical Dispersion Relation 


We now derive a linearized dispersion relation for the finite time step spectral algorithm, 
extending the arguments of §9-2 of [ 2 ] for our spectral algorithms. We begin by defining 
X = X -|- 5x, where X is the unperturbed ballistic motion, X^"") = x*^°) -|- and 

|k ■ 5x| <C 1 is a small perturbation. 

To include the boundary terms, we must include the update sequence from the 
end of the last “accelerate” to the beginning of the next. For convenience, we will 
denote these with integer indices, since they are all offset by a half-step. Inserting this 


approximation into eqn. (|12a|), along with our re-indexing, gives 

hei 




(n+l) 


+ 5x5"-'^ 


m 


/ dkkW(k)T(k-k.)(p, 


Jk-xF' 


(23) 


It is interesting to look at a finite width in T. This behaves like a finite width in v 
in the particle distribution, creating numerical Landau damping even for a cold plasma, 
with a damping time ~ ■ vq for a typical width Uk to the shape function T. We 

can intuitively understand this as follow: The finite width Fourier shapes represent a 
distribution of phase velocities. Thus, by introducing these shapes, we create a spread 
of phase velocities that a given mode couples to the distribution of particles, which is 
an unphysical but identical configuration to the usual mechanism for physical Landau 
damping. This may be useful in some scenarios. However, from here on we will consider 
only the discrete Fourier modes, \h(k — k„) = 5(k — kg-). 

Let us assume that the electrostatic field varies as (pa- = The integrand 

of the right hand side of this equation varies as ^nh)-iumh^ 

5x. Under this assumption, we have the dynamics for (5x given by 


and so too must 


5xf- 25x5") + 5x5"-') 


h 

m 


gikij ■ (x® nh) ^-iujnh 


(24) 
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It is convenient to define 


^x(") = 


which in turn gives that 


= --- - — 

2(cosa;2^h — 1) m 




(25) 


(26) 


It remains to compute $o-- Looking to eqn. (126) and Taylor expanding the right 
hand side gives 


\K\‘^K = 4:TTe^WjA*{K) 


Jkcy ■ nh) iumh- 


ikrr ■ SXj 


where we have assumed that the unperturbed quantity 

^TO^A*(k,)e‘''' ss 0 


(27) 


(28) 


which is to say that the unperturbed system has approximate charge neutrality. 

Combining the solution for the perturbed orbit with the Poisson equation yields 


~ ^ 

= dvr-^ WjA*{ko 


)e"*(k--ka')'X® X 


m 






2{cos u'^'h — 1) 

We note that the sum over j can be approximated as an integral over the initial 
conditions weighted by the local phase space density 




(30) 


where we have assumed the initial system has no charge perturbations. Noting the 
orthogonality relation with Fourier modes, we get that 


^° h^|A(k^)|^ [ dv/(v) 




(31) 


m J ' ' 2{cosuj^h — 1) 

where no is the unperturbed initial charge density. This yields the numerical dispersion 
relation 


1 = u}lh^\A{K)\^ / dv/(v)- 


4 sin 


2 ( {u}-'ka-v)h 


(32) 


for each mode k^-. In the limit of h —?• 0 and A = 1 (point particles) this returns the 
continuum limit. This is as in [ 2 ] but for discrete time steps but continuous space, 
although the discrete Fourier representation picks specific k^, and with the particle 
shape functions appearing to modify the dispersion relation. We can carry out additional 
analysis on this dispersion relation. 

Immediately we see that this algorithm modifies the plasma frequency for each 
mode to 


U}’= u}„\Mk,)\. 


(33) 
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This requires that |A(ko-)| stay close to unity until the Debye wavenumber, which is 
equivalent to requiring that the particle shape function A(x) be narrower than the 
Debye length. 

Assuming a cold coasting plasma, /(v) = (5(v — vq), the dispersion relation becomes 
4sin^ (34) 


which in turn requires that 

^|A(k,)| < 1 (35) 

which is the modihed Up — At condition for our electrostatic algorithm. This also gives a 
guideline for accurate spectral hdelity for the linear plasma dispersion relation (must be 
much smaller than 1), and the restriction is most stringent for long wavelength modes, 
assuming |A(ko-)| has a decreasing envelope with |ko-|. 

We can furthermore use Nyquist diagrams to evaluate the numerical Penrose 
criterion, since the dispersion relation has poles at a; — k^ ■ v = mr. However, only 
the n = 0 pole will contribute, as the residues of all the others are zero. Since dynamics 
is only parallel to k^., let us limit ourselves to a one-dimensional dispersion relation 
parallel to k^ for simplicity: 



This integral is given by the Cauchy principle value integral 

1 = ulh'^\A{K)\‘^ (v j dv ... + . (37) 

The imaginary part is of interest here. Because it has the identical functional form to 
the continuum limit with a modihed plasma frequency, we can conclude that there are 
no additional plasma instabilities introduced by the algorithm so long as the cUp — At 
condition is satished. 


That instability arises in the V part of the integral, which must be able to equal 
1 for real values of u. The sinusoid in the denominator may prevent that due to the 
repeated pole structure, although the details will depend on /(v). If we assume a cold 
but not inhnitely cold plasma, and a gaussian distribution f(v) = with 

V'2750-2 

(Jv V/ifc<T, then we can approximate the Cauchy principal value asymptotically as 



which means thermal effects offer a small modihcation of to the cold Up — At condition. 
It also offers an additional guideline for accurate modeling, which is that for an initially 


thermal distribution \hka<Jv\ 1 for the largest values of or else we could introduce 


a numerical plasma instability with a sufficiently large time step. 
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Figure 2: Total energy relative to the initial energy of the two-body system. 


6. Numerical Example — Two-Body Problem 

As a hrst test of the angular and linear momentum conserving properties of this 
algorithm, we consider the two dimensional Kepler problem. In two dimensions the 
particles are line charges, and the electrostatic potential is given by V{r) ~ ln(Y?’o), 
so this differs from the standard Yr Kepler problem. Our test problem is a single 
electron (mg, e) and a “double-positron” (2me,—2e). There is some initial total linear 
momentum. We use tent functions for the particle shapes and discrete Fourier modes, 
with T(k) = 5(k), with a time step of 5. x 10“®sec for 2 x 10® time steps. For the given 
initial conditions, this is approximately 9-10 steps per period, which reasonably resolves 
the dynamics of the problem. 

As is well-known, if one calculates the energy of a conservative Hamiltonian system 
from simulations using a symplectic integrator, that energy is bounded for exponential 
time. We see this behavior in £g. ([^. After some initial transient increase in the energy, 
it is conserved at the 1 % level over an extended time. We can also look at the individual 
energy contributions - kinetic, capacitive, and held energies - to see that each is very 
accurately bounded, so that there is no drift in the individual components that happens 
to cancel. The larger scale oscillations are a result of computing the energy every 250 
time steps, i.e. once every 25 periods or so. This creates much of the rippling we see in 

fig- (0- 

Fig. 0 demonstrates that the energy conservation does not arise from spurious 
growth in kinetic energy and decay in potential energy, or vice versa. As one 
would expect, the individual contributions oscillate and are well-behaved. The energy 
conservation then does not arise from any spurious trends in the individual contributions 
which happen to cancel out on the whole. 
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Figure 3: The individual contributions of the kinetic, capacitive, and held energies to 
the total energy of the Kepler problem. 
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Figure 4: Total Px and Py for the system, in a log scale. 


Furthermore, the robustness of the behavior is independent of the timing of the 
coordinates during the time step when the energy is computed. This means that the 
energy evaluated at the end of a full time step has the same level of precision (although 
will differ numerically) as the energy evaluated if a half-drift is used to put x and v 
at the same time, and then the helds are computed using x. This is another property 
familiar from single particle symplectic integrators. 

A major feature of this discrete Fourier spectral algorithm is that it is translationally 
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Figure 5: Energy error over a single time step versus h. 


invariant, therefore it should conserve the total momentum. In hg. (|^ we see that 
the initial total linear momentum is conserved to machine precision. There is an 
initial nonzero total momentum, and for the hrst hve hundred thousand time steps 
the momentum change is machine precision zero. The upward trend after this point is 
likely due to roundoff error. 

The algorithm also purports to be second order accurate, meaning that the 
correction to the continuous action integral is given hj Sjy ^ <S + 0{h^). Based on 
scaling, where S ^ f Edt, this would imply that the correction to the energy would be 
Ed k, E + 0{h‘^). For the two-body problem, we can clearly see this in £g. (|^. 

As can be seen, the energy deviation scales with for this particular problem 
from around 10“^ sec. where the failure to resolve the physics becomes catastrophicjjj] 
down to around 10“® sec where machine precision issues come into play. This clearly 
demonstrates the second order accuracy of the spectral algorithm. It should be noted 
that this algorithm was implemented using NumPy (v. 1.9.2) in Python (v. 2.7.10) 
on a 64-bit machine - handling of machine precision will vary across architectures and 
programming languages. 

We have thus shown that we can accurately simulate the two-body Kepler problem 
for a hundred thousand periods of the system with no change in the total momentum 
and with very well-behaved energy which does not spuriously grow. Such a problem is 
very difficult to do in conventional gridded PIC, where the total energy grows and the 
total momentum is changing, which breaks the careful conservation that allows for a 
closed orbit two-body problem. We see that this multisymplectic integrator preserves a 
number of key geometric properties of the system - the angular and linear momentum 

II At this point, the relative energy error is ^ 10^ — 10^. 







15 


Quantity 

Value [units] 

electron density 

10^^ cm“^ 

Ad 

3.96 X 10“^^ cm 

OJp 

8.14 X lO^^sec-i 

-^domain 

4 X Ad 

’^macro 

16 

macro, weight 

9.9 X 10® 

macro, width 

O.IAd 

6k 

V2(-^ciomain) 

^modes 

80 


5 X ‘^p/2-k 


Figure 6: First simulation parameters. 


are conserved exactly, with good behavior for the total energy due to the nature of 
symplectic algorithms. 

7. Numerical Example — Thermal Plasma 


Having demonstrated the energy and momentum conserving properties of this 
algorithm for a two-body problem, we now consider a thermal plasma. Again, we 
consider periodic boundaries to keep the system closed and make the total energy a 
quantity that is conserved in the physical system. In conventional PIC, such a system 
would show spurious energy growth over many plasma periods. 

We show simulation results run for ~ 10® plasma oscillations. We consider a system 
with the parameters in fig. (§, which is a very low resolution simulation. We are 
resolving Debye length physics here - in practice the spectral method introduces an 
ultraviolet cutoff on the resolved physics. Any length shorter than (fcmax)”^ is removed 
entirely, as can be verihed by failing to resolve the Debye length. 

In this simulation, there is no discernible change in the total energy of the system 
(fig- ( 0 )- This is to be expected, as the fields should be zero in the physical system, 
and we would expect all the motion to be ballistic. In two or three dimensions, we do 
not anticipate this level of energy conservation, which likely arises from a lack of degrees 
of freedom to cause the energy oscillations observed in the two-dimensional two-body 
problem. We also see the same momentum conservation, £g. (|^, which oscillates near 
machine precision throughout the simulation. 

In conventional particle-in-cell simulations, even in one dimension, it is known that 
the plasma will heat up [HI |T7] . As noted in na, various techniques in conventional 
algorithms - such as higher order particle shapes or increasing the grid resolution of the 
Debye length - can reduce the numerical heating up to a point. This algorithm removes 
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Figure 7: Energy change over time for thermal plasma simulation. 
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Figure 8: Momentum change over time for thermal plasma simulation. 


the heating entirely, as well as dramatically reducing the statistical fluctuations intrinsic 
in conventional particle-in-cell. 

The reduction in statistical fluctuations can be understood by looking at how the 
variational method would replicate a conventional algorithm. Conventional algorithms 
have no intrinsic convolution integrals, so we can understand the “particle shape” as 
either a particle shape with delta function helds, or a “held shape” with delta function 
particles. Both options are extremely noisy. By introducing the convolution integrals 
and using a Fourier basis, we smooth the particle shapes while simultaneously removing 
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high frequency content without aliasing. This has been demonstrated in finite-difference 
versions of variational algorithms in, e.g., [3], where a one dimensional electromagnetic 
finite differencing shows extremely smooth charge distributions for the applied problem. 
This comes from using a combination of field and particle shape functions, instead of 
just one or the other. 

8. Implementation and Optimization 


This algorithm has a number of properties that differ from conventional particle- 
in-cell methods that requires some consideration. This includes the discrete mode 
properties, parallel implementation, the global nature of the field solve, and evaluation 
of the deposition/interpolation terms. Also, the spectral nature of the algorithm requires 
some care to minimize the computational costs. 

The first is the span of fc-space. Implemented with periodic boundaries, the 
longest wavelength is L, the simulation domain length. Furthermore, we must resolve 
adequately the particle width, i. Thus, our /c-space modes must span k G [27r/L, ~ 
2n/i]. Furthermore, because the system has periodic boundaries, we must have as our 
increment in fc-space Ak = 27r/L to maintain periodicity. Thus, we must use 0{L/i)^ 
Fourier modes. This is comparable to the number of cells required for a traditional 
gridded method. 

Because the algorithm is truly spectral (there is no grid in real space) conventional 
domain decomposition is irrelevant. This simplifies parallel implementations. Each 
process must have all the field data, and therefore the field solve can be reduced to a 
single MPI_AllReduce call on the after the deposition phase. All other operations can 
be performed locally. This also makes load balancing trivial - every process will receive 
an equal number of macro-particles and store all the field information. However, this 
could present a challenge for very large simulations with substantial amounts of field 
data. Because the simulation requires no overhead for messaging particle data, and is 
perfectly load balanced, such spectral algorithms may be ideal for simulations of devices 
such as magnetrons which have fields everywhere but particles which are isolated to a 
small region of space throughout much of the simulation. 

Computing the deposition and interpolation can be computationally intensive due 
to the exponential evaluations. This can be optimized in two ways for the discrete 
Fourier treatment (Fourier shape functions cannot avoid this problem). The first is by 
recognizing that we are always taking the real part of the fields, so we only have to keep 
the k„ with all components greater than or equal to zero. This is a factor of 2^ savings 
in the amount of information that must be stored. We can also reduce the number of 
exponentials per particle from Nmodes to D if we use regularly spaced Fourier modes 
spaced by by tracking Then we must evaluate the exponential once, and each 

coefficient is simply a positive integer power of that exponential. Compared to the the 
naive algorithm, this can represent a substantial savings in computation and memory 
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footprint. 

9. Conclusion 

We have derived an electrostatic spectral macroparticle algorithm from a variational 
approach. The resulting algorithm closely resembles conventional leapfrog algorithms, 
although we have bypassed using a grid - hence the choice to refer to it as a “spectral 
macroparticle algorithm” instead of a “particle in cell algorithm”. A broad class of 
algorithms are available with hnite widths in /c-space for the field modes. However, 
our analysis has shown that discrete Fourier modes (rather than modes with hnite 
width) are exactly momentum conserving, and furthermore observed that hnite width 
modes introduce a numerical Landau damping term to the dispersion relation. For these 
reasons, we consider the discrete Fourier algorithm as the ideal implementation. This 
algorithm has demonstrated remarkable energy conservation properties, as one would 
expect from symplectic algorithms. 

We have demonstrated the discrete Fourier algorithm prevents all spatial 
discretization instabilities, and is only prone to the Up — At instability that arises 
from time discretization. This is therefore a remedy for the hnite grid instability. We 
have also seen these algorithms demonstrated in two two-dimensional simulations: a 
two-body problem and a thermal plasma, both with periodic boundaries. These are, 
to the author’s best knowledge, the hrst computational examples of two-dimensional 
variational algorithms. They demonstrate proper second-order energy error, as one 
would expect. 

Because we considered Cartesian electrostatics, we were not presented with a 
number of issues important for electromagnetic algorithms, particularly gauge invariance 
and the associated charge conservation, as well as difficulties in non-Cartesian coordinate 
systems such as the pole at the origin of spherical or cylindrical coordinates. This 
represents future work. 

Variational algorithms have a number of advantages over conventional plasma 
algorithms. The primary one is that the equations of motion come from a discrete 
Euler-Lagrange equation. Because of this, all approximations (spectral decomposition, 
coordinate transformations for particles and helds, time discretization) are made at the 
level of the action - the discrete Euler-Lagrange equations then construct the algorithm 
from these approximations. This gives more freedom to the algorithm developer, as the 
approximations are guaranteed to be self-consistent. 
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